rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra,readxl)
setwd('put_your_wd_here')


# shootings
allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date <- mdy(allshootings$date)
top10 <- allshootings$shooting[allshootings$top10==1]

filloutdates <- function(x){
  x <- left_join(
    data.frame(donationdate = seq(range(x$donationdate)[1],range(x$donationdate)[2], by='day')),
    x, by='donationdate'
  ) 
  x$donationcount[is.na(x$donationcount)] <- 0
  x$contribution_receipt_amount[is.na(x$contribution_receipt_amount)] <- 0
  x$Group <- x$Group[1]
  return(x)
}

giffords <- read_stata('datasets/GiffordsPACDailyDonations.dta')
nra <- read_stata('datasets/NRADailyDonations.dta')

giffords <- filloutdates(giffords)
nra <- filloutdates(nra)

filtDates <- function(df){
  df$date <- ymd(df$donationdate)
  mindate <- min(df$donationdate) + months(2)
  maxdate <- max(df$donationdate) - months(2)
  dates <- allshootings %>%
    filter(date >= mindate & date <= maxdate & include==1)
  return(dates)
}

bottom20 <- read_csv('datasets/lower_20_donations.csv')
  
# RDDTOOLS

makeEstRDD <- function(out, treat, name, org, time=2){
  
  min <- treat - months(time) 
  max <- treat + months(time) 
  out$running_day <- as.numeric(out$donationdate - treat)
  
  sd_est <- sd(out$donationcount,na.rm=T)
  house_rdd <- rdd_data(y=out$donationcount, x=out$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name,
                        outcome=org)
  return(fullest)
}

giffordsdates <- filtDates(giffords)
nradates <- filtDates(nra)

gifford_results <- vector('list', length(giffordsdates$date))
for(i in 1:length(gifford_results)){
  gifford_results[[i]] <- makeEstRDD(giffords, giffordsdates$date[i], name = giffordsdates$shooting[i], org='Giffords PAC')
  print(i)
}

nra_results <- vector('list', length(nradates$date))
for(i in 1:length(nra_results)){
  nra_results[[i]] <- makeEstRDD(nra, nradates$date[i], name = nradates$shooting[i], org='NRA PAC')
  print(i)
}

dat <- bind_rows(
  bind_rows(gifford_results),
  bind_rows(nra_results)
) %>%
  mutate()

dat$group_var <- 'Number Donations'
dat1 <- dat

makeEstAmntRDD <- function(out, treat, name, org,time=2){
  min <- treat - months(time)
  max <- treat + months(time)
  out <- filter(out, 
                donationdate >= min & donationdate <= max)
  out$running_day <- as.numeric(out$donationdate - treat)
  out <<- out
  
  sd_est <- sd(out$contribution_receipt_amount, na.rm=T)
  house_rdd <- rdd_data(y=out$contribution_receipt_amount, x=out$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name,
                        outcome=org)
  return(fullest)
}

gifford_results <- vector('list', length(giffordsdates$date))
for(i in 1:length(gifford_results)){
  gifford_results[[i]] <- makeEstAmntRDD(giffords, giffordsdates$date[i], name = giffordsdates$shooting[i], org='Giffords PAC')
  print(i)
}

nra_results <- vector('list', length(nradates$date))
for(i in 1:length(nra_results)){
  nra_results[[i]] <- makeEstAmntRDD(nra, nradates$date[i], name = nradates$shooting[i], org='NRA PAC')
  print(i)
}

dat <- bind_rows(
  bind_rows(gifford_results),
  bind_rows(nra_results)
)

dat$group_var <- 'Donations Amount'

dat2 <- dat

uniqshoots <- allshootings$shooting[allshootings$include==1] %>% na.omit()
empty <- expand.grid(uniqshoots, c('Donations Amount','Number Donations'), c('Giffords PAC','NRA PAC')) %>%
  mutate(Coef=NA, lower=NA, upper=NA) %>%
  rename(shooting = Var1, group_var = Var2, outcome=Var3) %>%
  filter(!shooting %in% dat1$shooting)

dat <- bind_rows(dat1,dat2, empty) 

write_csv(dat, file='datasets/donations_meta.csv')

####

allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date <- mdy(allshootings$date)

dat <- read_csv('datasets/donations_meta.csv')

dat$statsig <- ifelse(dat$lower > 0 | dat$upper < 0,1,0)
dat$statsig[is.na(dat$statsig)] <- 0
dat$statsig <- factor(dat$statsig)
dat$shooting <- factor(dat$shooting, levels=rev(c(
  allshootings$shooting
)))

calcMeta <- function(data, name, outcome, shooting, group_var){
  m.gen <- meta::metagen(TE = Coef,
                   seTE = se,
                   studlab = shooting,
                   data = data,
                   sm = "SMD",
                   fixed = FALSE,
                   random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    outcome = outcome,
                    group_var = group_var,
                    statsig=factor(1))
  return(results)
}

datsub <- dat %>% 
  filter(outcome %in% c('Giffords PAC'))

meta_number <- calcMeta(datsub %>% filter(group_var == 'Number Donations') %>% mutate(se = (upper - Coef)/1.96), 
                 'Meta Estimate (Top 20)','Giffords PAC','Meta Estimate (All)', 'Meta Number') %>%
  mutate(statsig=factor(0))
meta_amount <- calcMeta(datsub %>% filter(group_var == 'Donations Amount') %>% mutate(se = (upper - Coef)/1.96), 
                        'Meta Estimate (Top 20)','Giffords PAC','Meta Estimate (All)', 'Meta Amount') %>%
  mutate(statsig=factor(1))
meta_number2 <- calcMeta(datsub %>% filter(group_var == 'Number Donations' & shooting %in% top10) %>% mutate(se = (upper - Coef)/1.96), 
                        'Meta Estimate (Top 10)','Giffords PAC','Meta Estimate (Top 10)', 'Meta Number') %>%
  mutate(statsig=factor(0))
meta_amount2 <- calcMeta(datsub %>% filter(group_var == 'Donations Amount' & shooting %in% top10) %>% mutate(se = (upper - Coef)/1.96), 
                        'Meta Estimate (Top 10)','Giffords PAC','Meta Estimate (Top 10)', 'Meta Amount') %>%
  mutate(statsig=factor(1))
bottom20$shooting <- 'Meta Estimate (Bottom 20)'
bottom20 <- bottom20 %>% mutate(statsig=factor(0))

bind_rows(meta_amount,meta_number,
          meta_amount2, meta_number2, bottom20[1:2,]) %>%
  mutate(z='Donation Giffords Meta Analyses') %>%
  write_csv(file='datasets/meta_donations_giffords.csv')


p1 <- bind_rows(datsub,meta_amount,meta_number,
                meta_amount2, meta_number2, bottom20[1:2,]) %>%
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting,'Meta Estimate (Top 10)',
    'Meta Estimate (Top 20)', 'Meta Estimate (Bottom 20)'
  )))) %>%
  mutate(group_var=factor(group_var,levels=c(
    'Number Donations','Donations Amount',
    'Meta Number','Meta Amount'
  ))) %>%
  ggplot(aes(shooting, y=Coef, ymin=lower, ymax=upper, 
           color=statsig, shape=group_var)) +
    coord_flip() +
    geom_hline(yintercept=0, linetype=2, color='red') +
    labs(title='(A) Gun Control\nDonations', x='',y='') +
    theme_bw() + 
    geom_errorbar(width=0, position=position_dodge(width=0.5)) +
    geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
    scale_color_manual(values=c('grey50','Black')) +
    scale_shape_manual(values = c(24, 21,17,16)) + 
  facet_wrap(~outcome) +
    theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5)) 
  
p1

datsub <- dat %>% filter(outcome %in% c('NRA PAC'))
meta1 <- calcMeta(datsub %>% filter(group_var == 'Number Donations') %>% mutate(se = (upper - Coef)/1.96), 
                   'Meta Estimate (Top 20)','NRA PAC','Meta Analysis', 'Meta Number')%>%
  mutate(statsig=factor(0))
meta2 <- calcMeta(datsub %>% filter(group_var == 'Donations Amount') %>% mutate(se = (upper - Coef)/1.96), 
                 'Meta Estimate (Top 20)','NRA PAC','Meta Analysis', 'Meta Amount') %>%
  mutate(statsig=factor(0))
meta3 <- calcMeta(datsub %>% filter(group_var == 'Number Donations' & shooting %in% top10) %>% mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (Top 10)','NRA PAC','Meta Analysis', 'Meta Number')%>%
  mutate(statsig=factor(0))
meta4 <- calcMeta(datsub %>% filter(group_var == 'Donations Amount' & shooting %in% top10) %>% mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (Top 10)','NRA PAC','Meta Analysis', 'Meta Amount') %>%
  mutate(statsig=factor(0))

bind_rows(meta1, meta2, meta3, meta4, bottom20[3:4,]) %>%
  mutate(z='Donation NRA Meta Analyses') %>%
  write_csv(file='datasets/meta_donations_NRA.csv')

p2 <- bind_rows(datsub,meta1, meta2, meta3, meta4, bottom20[3:4,]) %>%
  mutate(group_var = case_when(
    group_var == 'Number Donations' ~ 'Number',
    group_var == 'Donations Amount' ~ 'Amount',
    TRUE ~ group_var
  ),
    group_var = factor(group_var, levels=c(
    'Number','Amount', 'Meta Number', 'Meta Amount'
  ))) %>%
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting,'Meta Estimate (Top 10)','Meta Estimate (Top 20)',
    'Meta Estimate (Bottom 20)'
  )))) %>%
  ggplot(aes(shooting, y=Coef, ymin=lower, ymax=upper, color=statsig, shape=group_var)) +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(B) Gun Rights\nDonations', x='',y='') +
  theme_bw() + 
  geom_errorbar(width=0, position=position_dodge(width=0.5)) +
  geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
  scale_color_manual(values=c('grey50','Black'), guide='none') + 
  scale_shape_manual(values = c(21, 24, 16, 17), name='') + 
  facet_wrap(~outcome) +
  theme(plot.title = element_text(hjust = 0.5), 
        axis.text.y = element_blank(),
        legend.position = 'none') 
p2

plots <- ggarrange(p1,p2,
                   heights = c(4,4),
                   widths = c(6,4),
                   ncol = 2, nrow = 1)
plots
plots <- annotate_figure(plots,
                bottom = text_grob("Effect on Donations\n(Standard Deviations)",family='Helvetica')) 
ggsave(plots, width=7, height=5, filename = 'donations_rddtools.png')


####
# Meta analyses by race of victims
#### 

race_est <- read_csv('datasets/shooting_race.csv')
race_est <- race_est %>%
  filter(shooting %in% dat$shooting)
terciles <- quantile(race_est$mean.white, na.rm=T, c(0,0.33,0.66,1))
race_est$tercile <- ifelse(race_est$mean.white < terciles[2],1,0)
race_est$tercile <- ifelse(race_est$mean.white > terciles[2] & race_est$mean.white < terciles[3],2,race_est$tercile)
race_est$tercile <- ifelse(race_est$mean.white > terciles[3],3,race_est$tercile)
datnew <- left_join(dat, race_est, by='shooting')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & outcome == 'NRA PAC' & group_var == 'Number Donations'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis','x')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & outcome == 'NRA PAC' & group_var == 'Number Donations'),  
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis','x')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & outcome == 'NRA PAC' & group_var == 'Number Donations'),  
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis','x')
meta_tercile_nra_number <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Donations', z = 'NRA', y = 'Number')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & outcome == 'NRA PAC' & group_var == 'Donations Amount'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis','x')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & outcome == 'NRA PAC' & group_var == 'Donations Amount'),  
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis','x')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & outcome == 'NRA PAC' & group_var == 'Donations Amount'),  
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis','x')
meta_tercile_nra_amount <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Donations', z = 'NRA', y = 'Amount')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & outcome == 'Giffords PAC' & group_var == 'Number Donations'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis','x')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & outcome == 'Giffords PAC' & group_var == 'Number Donations'),  
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis','x')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & outcome == 'Giffords PAC' & group_var == 'Number Donations'),  
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis','x')
meta_tercile_giffords_number <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Donations', z = 'Giffords', y = 'Number')

meta_top3 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 3 & outcome == 'Giffords PAC' & group_var == 'Donations Amount'),  
                      'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis','x')
meta_top2 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 2 & outcome == 'Giffords PAC' & group_var == 'Donations Amount'),  
                      'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis','x')
meta_top1 <- calcMeta(datnew %>% 
                        mutate(se = (upper - Coef)/1.96) %>% 
                        filter(tercile == 1 & outcome == 'Giffords PAC' & group_var == 'Donations Amount'),  
                      'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis','x')
meta_tercile_giffords_amount <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
  mutate(category = 'Donations', z = 'Giffords', y = 'Amount')


dat.out <- bind_rows(
  meta_tercile_nra_number, 
  meta_tercile_nra_amount,
  meta_tercile_giffords_number,
  meta_tercile_giffords_amount
)
write_csv(dat.out, 'datasets/metaanalysis/meta_race_donations.csv' )


# Bottom 20 -----
rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra,readxl)
setwd('put_your_wd_here')

# shootings
allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date <- mdy(allshootings$date)

filloutdates <- function(x){
  x <- left_join(
    data.frame(donationdate = seq(range(x$donationdate)[1],range(x$donationdate)[2], by='day')),
    x, by='donationdate'
  ) 
  x$donationcount[is.na(x$donationcount)] <- 0
  x$contribution_receipt_amount[is.na(x$contribution_receipt_amount)] <- 0
  x$Group <- x$Group[1]
  return(x)
}
giffords <- read_stata('datasets/GiffordsPACDailyDonations.dta')
nra <- read_stata('datasets/NRADailyDonations.dta')

giffords <- filloutdates(giffords)
nra <- filloutdates(nra)

filtDates <- function(df){
  df$date <- ymd(df$donationdate)
  mindate <- min(df$donationdate) + months(2)
  maxdate <- max(df$donationdate) - months(2)
  dates <- allshootings %>%
    filter(date >= mindate & date <= maxdate & include==0)
  return(dates)
}

makeEstRDD <- function(out, treat, name, org, time=2){
  
  min <- treat - months(time) 
  max <- treat + months(time) 
  out$running_day <- as.numeric(out$donationdate - treat)
  
  sd_est <- sd(out$donationcount,na.rm=T)
  house_rdd <- rdd_data(y=out$donationcount, x=out$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name,
                        outcome=org)
  return(fullest)
}

giffordsdates <- filtDates(giffords)
nradates <- filtDates(nra)

gifford_results <- vector('list', length(giffordsdates$date))
for(i in 1:length(gifford_results)){
  gifford_results[[i]] <- makeEstRDD(giffords, giffordsdates$date[i], name = giffordsdates$shooting[i], org='Giffords PAC')
  print(i)
}

nra_results <- vector('list', length(nradates$date))
for(i in 1:length(nra_results)){
  nra_results[[i]] <- makeEstRDD(nra, nradates$date[i], name = nradates$shooting[i], org='NRA PAC')
  print(i)
}

dat <- bind_rows(
  bind_rows(gifford_results),
  bind_rows(nra_results)
) %>%
  mutate()

dat$group_var <- 'Number Donations'
dat1 <- dat

makeEstAmntRDD <- function(out, treat, name, org,time=2){
  min <- treat - months(time)
  max <- treat + months(time)
  out <- filter(out, 
                donationdate >= min & donationdate <= max)
  out$running_day <- as.numeric(out$donationdate - treat)
  out <<- out
  
  sd_est <- sd(out$contribution_receipt_amount, na.rm=T)
  house_rdd <- rdd_data(y=out$contribution_receipt_amount, x=out$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name,
                        outcome=org)
  return(fullest)
}


gifford_results <- vector('list', length(giffordsdates$date))
for(i in 1:length(gifford_results)){
  gifford_results[[i]] <- makeEstAmntRDD(giffords, giffordsdates$date[i], name = giffordsdates$shooting[i], org='Giffords PAC')
  print(i)
}

nra_results <- vector('list', length(nradates$date))
for(i in 1:length(nra_results)){
  nra_results[[i]] <- makeEstAmntRDD(nra, nradates$date[i], name = nradates$shooting[i], org='NRA PAC')
  print(i)
}

dat <- bind_rows(
  bind_rows(gifford_results),
  bind_rows(nra_results)
)

dat$group_var <- 'Donations Amount'
dat2 <- dat

uniqshoots <- allshootings$shooting[allshootings$include==0] %>% na.omit()
empty <- expand.grid(uniqshoots, c('Donations Amount','Number Donations'), c('Giffords PAC','NRA PAC')) %>%
  mutate(Coef=NA, lower=NA, upper=NA) %>%
  rename(shooting = Var1, group_var = Var2, outcome=Var3) %>%
  filter(!shooting %in% dat1$shooting)

dat <- bind_rows(dat1,dat2, empty) 

#write_csv(dat, file='datasets/donations_meta.csv')

dat$statsig <- ifelse((dat$lower > 0 & dat$upper > 0) | (dat$lower < 0 & dat$upper <0),1,0)
dat$statsig[is.na(dat$statsig)] <- 0
dat$statsig <- factor(dat$statsig)
dat$shooting <- factor(dat$shooting, levels=rev(c(
  allshootings$shooting
)))

calcMeta <- function(data, name, outcome, shooting, group_var){
  m.gen <- meta::metagen(TE = Coef,
                         seTE = se,
                         studlab = shooting,
                         data = data,
                         sm = "SMD",
                         fixed = FALSE,
                         random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    outcome = outcome,
                    group_var = group_var,
                    statsig=factor(1))
  return(results)
}

datsub <- dat %>% 
  filter(outcome %in% c('Giffords PAC'))

meta_number <- calcMeta(datsub %>% filter(group_var == 'Number Donations') %>% mutate(se = (upper - Coef)/1.96), 
                        'Meta Estimate (All)','Giffords PAC','Meta Estimate (All)', 'Meta Number') %>%
  mutate(statsig=factor(0))
meta_amount <- calcMeta(datsub %>% filter(group_var == 'Donations Amount') %>% mutate(se = (upper - Coef)/1.96), 
                        'Meta Estimate (All)','Giffords PAC','Meta Estimate (All)', 'Meta Amount') %>%
  mutate(statsig=factor(0))

p1 <- bind_rows(datsub,meta_amount,meta_number) %>%
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting,'Meta Estimate (All)'
  )))) %>%
  ggplot(aes(shooting, y=Coef, ymin=lower, ymax=upper, color=statsig, shape=group_var)) +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(B) Gun Control\nDonations', x='',y='') +
  theme_bw() + 
  geom_errorbar(width=0, position=position_dodge(width=0.5)) +
  geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
  theme(legend.position='none', 
        plot.title = element_text(hjust = 0.5)) +
  scale_color_manual(values=c('grey50','Black')) +
  scale_shape_manual(values = c(24, 17, 16, 21)) + 
  facet_wrap(~outcome)
p1

datsub <- dat %>% filter(outcome %in% c('NRA PAC'))
meta1 <- calcMeta(datsub %>% filter(group_var == 'Number Donations') %>% mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (All)','NRA PAC','Meta Analysis', 'Meta Number')%>%
  mutate(statsig=factor(0))
meta2 <- calcMeta(datsub %>% filter(group_var == 'Donations Amount') %>% mutate(se = (upper - Coef)/1.96), 
                  'Meta Estimate (All)','NRA PAC','Meta Analysis', 'Meta Amount') %>%
  mutate(statsig=factor(0))

# write out metas 
bind_rows(meta_number, meta_amount, meta1, meta2) %>%
  write_csv(file = 'datasets/lower_20_donations.csv')

p2 <- bind_rows(datsub,meta1, meta2) %>%
  mutate(group_var = case_when(
    group_var == 'Number Donations' ~ 'Number',
    group_var == 'Donations Amount' ~ 'Amount',
    TRUE ~ group_var
  ),
  group_var = factor(group_var, levels=c(
    'Number','Amount', 'Meta Number', 'Meta Amount'
  ))) %>%
  mutate(shooting = factor(shooting, levels=rev(c(
    allshootings$shooting,'Meta Estimate (All)'
  )))) %>%
  ggplot(aes(shooting, y=Coef, ymin=lower, ymax=upper, color=statsig, shape=group_var)) +
  coord_flip() +
  geom_hline(yintercept=0, linetype=2, color='red') +
  labs(title='(C) Gun Rights\nDonations', x='',y='') +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5), 
        axis.text.y = element_blank(),
        legend.position = 'none') +
  geom_errorbar(width=0, position=position_dodge(width=0.5)) +
  geom_point(position=position_dodge(width=0.5), size=2, fill='white') +
  scale_color_manual(values=c('grey50','Black'), guide='none') + 
  scale_shape_manual(values = c(21, 24, 16, 17), name='') + 
  facet_wrap(~outcome)
p2

plots <- ggarrange(p1,p2,
                   heights = c(4,4),
                   widths = c(6,4),
                   ncol = 2, nrow = 1)
plots

ggsave(plots, width=7, height=5, filename = 'figures/donations_rddtools_lowatten.png')


